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Popular Summary 

Global maps of rainfall are of great importance in connection with modeling of the earth’s climate. 
Comparison between the maps of rainfall predicted by computer-generated climate models with 
observation provides a sensitive test for these models. To make such a comparison, one typically needs 
the total precipitation amount over a large area, which could be hundreds of kilometers in size over 
extended periods of time of order days or months. This presents a difficult problem since rain varies 
greatly from place to place as well as in time. 

Remote sensing methods using ground radar or satellites detect rain over a large area by essentially taking 
a series of snapshots at infrequent intervals and indirectly deriving the average rain intensity within a 
collection of “pixels”, usually several kilometers in size. They measure area average of rain at a particular 
instant. Rain gauges, on the other hand, record rain accumulation continuously in time but only over a 
very small area tens of centimeters across, say, the size of a dinner plate. They measure only a time 
average at a single location. In making use of either method one needs to fill in the gaps in the observation 
- either the gaps in the area covered or the gaps in time of observation. This involves using statistical 
models to obtain information about the rain that is missed from what is actually detected. This paper 
investigates such a statistical model and validates it with rain data collected over the tropical Western 
Pacific from ship bome radars during TOGA CO ARE (Tropical Oceans Global Atmosphere Coupled 
Ocean-Atmosphere Response Experiment). 

The model incorporates a number of commonly observed features of rain. While rain varies rapidly with 
location and time, the variability diminishes when averaged over larger areas or longer periods of time. 
Moreover, rain is patchy in nature - at any instant on the average only a certain fraction of the observed 
pixels contain rain. The fraction of area covered by rain decreases, as the size of a pixel becomes smaller. 
This means that within what looks like a patch of rainy area in a coarse resolution view with larger pixel 
size, one finds clusters of rainy and dry patches when viewed on a finer scale. The model makes definite 
predictions about how these and other related statistics depend on the pixel size. These predictions were 
found to agree well with data. In a subsequent second part of the work we plan to test the model with rain 
gauge data collected during the TRMM (Tropical Rainfall Measuring Mission) ground validation 
campaign. 
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Abstract. A characteristic feature of rainfall statistics is that they depend 
on the space and time scales over which rain data are averaged. As part of 
an earlier effort to determine the sampling error of satellite rain averages, a 
space-time model of rain statistics was developed to describe the statistics of 
gridded rain observed over the eastern tropical Atlantic. The model allows 
one to compute the second moment statistics of space- and time-averaged rain 
rate, which can be fitted to satellite or rain gauge data to determine the four 
model parameters - an overall strength parameter, a characteristic length 
separating the long and short wavelength regimes, a characteristic relaxation 
time for decay of the autocorrelation of the instantaneous local rain rate, and 
a certain “fractal” power law exponent. For area-averaged instantaneous rain 
rate, this exponent governs the power law dependence of these statistics on 
the averaging length scale L. In particular, the variance of rain rate averaged 
over an L x L area exhibits a power law singularity as L — » 0. In the present 
work a more efficient method of estimating the model parameters is presented, 
and used to fit the model to the statistics of area-averaged rain rate over 
the tropical Western Pacific measured with ship-borne radar during TOGA 
COARE. Good agreement is found between the data and predictions from 
the model over a wide range of averaging scales. An extension of the spectral 
model scaling relations to describe the dependence of the fractional coverage 
of rain on the spatial scale L is also explored. 

1. Introduction 

Rainfall is a complex phenomenon involving the in- 
terplay of many physical processes in the atmosphere. 

This makes it extremely difficult to quantitatively de- 
scribe rainfall in terms of a detailed model of the un- 
derlying processes. Instead, a mathematically tractable 
and yet realistic description naturally calls for statisti- 
cal considerations. In contrast with the physical mod- 
els, models of rain statistics have the advantage of be- 
ing conceptually economical in the sense that they gen- 
erally involve fewer adjustable parameters. Moreover, 


they can be easily validated by comparing against real 
data gathered over a large area and a large time pe- 
riod representative of a certain climatological regime. 
Once the model parameters are determined from a suf- 
ficiently large data set, the model provides an efficient 
method of describing various statistical properties of 
rainfall over areas with similar rain climatologies [see, 
e.g., Bell et al, 2001, hereafter referred to as BKKj. For 
many purposes, this is all that is needed. 

It is common experience that rainfall is highly ir- 
regular in space and time and difficult to accurately 
predict and measure. Rain gauges have traditionally 
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provided a direct measurement of rain on a local ba- 
sis. However, the inherently sporadic nature of rainfall 
along with the practical difficulty of laying down a suffi- 
ciently dense network of continuously monitored gauges 
makes this method ill-suited to the problem of measur- 
ing precipitation on a global scale over an extended time 
period. In recent years, remote sensing techniques have 
provided an alternative and more convenient method 
of rain estimation. These techniques are attractive be- 
cause of their ability to cover large areas at a time. Re- 
mote sensing measurements, when carried out from low 
earth-orbiting satellites, enable one to estimate rainfall 
over a truly global scale. 

Inherent in these remote sensing measurements is a 
finite spatial resolution, since these methods infer the 
rain rate indirectly from the intensity of the radiation 
that carries information about precipitation present 
within a certain volume of the atmosphere. This resolu- 
tion is governed by the spatial resolution of the instru- 
ment’s antennae. Thus, the primary quantity estimated 
by these methods from the retrieval output is the area- 
average of the local instantaneous rain rate i?(x, t) over 
a spatial domain A of area A , which approximates the 
effective area covered by the instrument (whose actual 
response profile may be quite complex), 

K A (t) = jfj 2 x R (xi,t). (1) 

It is, however, somewhat of an abstraction, since the ac- 
tual radiative information resulting from the measure- 
ment process is more closely related to an area- and 
time-average over a finite duration of the order of the 
time the rain takes to fall to the ground. By contrast, a 
rain gauge provides a measurement of the time average 
at a point x (to the extent a rain gauge can be treated 
as a geometrical point), over a finite time T: 

i r T/2 

72r(x,$) = ^ / drRfct + r). ( 2 ) 

1 J-T/2 

It is worth noting that both techniques of measurement 
are affected by systematic as well as random errors. 
While remote sensing estimates generally suffer from 
retrieval errors due to uncertainties in the rain retrieval 
algorithms used to convert the radiances into rain rates, 
the rain gauges provide a direct measure of the amount 
of rain at a location. However, rain gauge data can con- 
tain gaps and may be unreliable due to difficulties in 
continuous monitoring and occasional instrument fail- 
ure. In addition, rain gauge data are rather sensitive to 
local factors, such as wind and evaporation, orographic 
details and placement of the gauges. 


Rainfall, despite its apparently irregular nature, is 
correlated both in space and in time. Spatio-temporal 
variability of rain directly influences other hydrologi- 
cal processes, such as soil moisture transport and sur- 
face run-off. Statistical behavior of area- and/or time- 
averaged rain rate at different length and time scales 
can also be an important diagnostic for comparing ac- 
tual rain behavior with predictions from global climate 
models (which in general only predict area-averaged 
rain rate over relatively large grid boxes, typically any- 
where from 10 km to several degrees). They also have 
important consequences for the estimation of the so- 
called sampling error that arises in satellite measure- 
ments because of the intermittent and often incomplete 
coverage of a given area of the globe [McConnell and 
North , 1987; Bell et a/., 1990; Bell and Kundu , 1996, 
hereafter BK96]. 

In this paper we present a statistical model of space 
time variability of tropical rainfall seen in precipitation 
data sets created through ground based radar observa- 
tions. We describe the data in terms of a simple model 
used previously by BK96 in their study of the sampling 
error in preparation for the Tropical Rainfall Measure- 
ment Mission (TRMM). In this model Fourier modes of 
the local instantaneous rain rate i?(x, t) are treated as a 
set of random variables whose evolution is governed by 
first order linear stochastic equations driven by white 
noise sources. This equation immediately leads to an 
analytical form of the rainfall power spectrum depend- 
ing on a few adjustable parameters. The Fourier trans- 
form of the spectrum yields the space-time covariance 
function of the local rain rate field, 

c(x,t;x',i') = (R'(x,t)R'(x',t')) , (3) 

where 

i?'(x, t) = i?(x, t) - (i?(x, t)) (4) 

and the angle brackets represent a statistical average 
over a hypothetical collection, or ensemble of rain fields 
with similar rain climatology. This function when in- 
tegrated over finite grid box areas or over finite time 
intervals with respect to each of its space or time argu- 
ments, yields the space-time covariance function of the 
area- or time-averaged rain rate, which are then fitted 
to the data to determine the model parameters. The 
model can easily address questions with regard to de- 
pendence of the statistical properties of rain on spatial 
and temporal resolution scales. This model, which in 
BK96 was fairly successful in describing the Global At- 
mospheric Research Program (GARP) Atlantic Trop- 
ical Experiment (GATE) Phase I data, is here tested 
with a gridded precipitation data set constructed from 
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the radar images obtained during the Tropical Ocean 
Global Atmosphere (TOGA) Coupled Ocean Atmo- 
sphere Response Experiment (COARE) [Webster and 
Lukas , 1992]. This data set covers a longer period of 
observation than the GATE data set studied in BK96 
and is gathered from a region in the tropical western Pa- 
cific rather than the eastern Atlantic where GATE was 
conducted. Since the data set gives us instantaneous 
spatial averages on a regular rectangular cartesian grid 
analogous to the GATE data, it is suitable for studying 
the statistical behavior of the rain field averaged over 
various spatial scales. 

There exists a large class of physical rainfall models 
based on random processes that may be called hier- 
archical cluster models. These models attempt to cap- 
ture the spatial and temporal organization of storm sys- 
tems that cause precipitation. A general mathematical 
framework for such a model was laid by Le Cam [1961] 
in terms of “marked point processes” . He used stochas- 
tic Poisson processes with suitably randomized param- 
eters to describe the irregular occurrence of rain events 
in space and time. In many models spatial aggrega- 
tion is described in terms of an organized hierarchy of 
structures - “rainbands” of elongated elliptical shape 
with circular rain “cells” clustered within them. The 
“marks” refer to a set of physical attributes which may 
include maximum rainfall intensity, size and speed of 
the cells and so on. A number of such models of var- 
ious degrees of complexity exist in the literature [e.g., 
Gupta and Waymire , 1979; Cox and Isham , 1988; Smith 
and Krajewski , 1987]. A particularly elaborate stochas- 
tic space-time model of precipitation with exponentially 
decaying gaussian distributed rain intensity field within 
a rain cell was introduced by Waymire et al [1984] fol- 
lowing earlier work by Gupta and Waymire [1979]. The 
Waymire-Gupta-Rodrfguez-Iturbe (WGR) model [see 
Waymire et al , 1984] aims to describe the mesoscale 
structure of precipitation in some detail, including for 
instance the effect of global advection of rain bands, at 
the price of introducing a large number of phenomeno- 
logical parameters. Estimation of the parameters char- 
acterizing a cluster model from actual rain data poses 
a difficult problem. Parameter estimation for simpler 
cluster models has been investigated by Smith and Karr 
[1987], Valdes et al. [1985], Smith and Krajewski [1987] 
and Krajewski and Smith [1989]. Islam et al [1988], 
Valdes et al [1990], and Koepsell and Valdes [1991] used 
nonlinear least squares fitting to estimate the WGR 
model parameters - a difficult task in view of the com- 
plexities of the model. An analytical form of the model 
spectrum was derived by Valdes et al [1990]. However, 


Valdes et al [1994] have found that the WGR model 
spectrum using the estimated model parameters does 
not describe the GATE spectrum well. For a recent ac- 
count of a number of models in this category, including 
the problem of parameter estimation in the context of a 
precipitation data set obtained in the United Kingdom 
from the HYREX project, a large network of radar and 
rain gauges, see Wheater et al [2000]. It is the ease of 
parameter estimation and simplicity of application to 
real rain data sets, combined with its ability to capture 
the variability on many different space and time scales, 
that motivated us to consider a rainfall model based on 
a simple stochastic dynamical equation rather than a 
cluster model of the type described above. 

Other rainfall models based on a stochastic dy- 
namical equation exist in the literature [North and 
Nakamoto , 1989; Yoo et al , 1996] which also describe 
the GATE spectrum with varying degrees of success. Of 
these the North-N akamoto (NN) forced diffusion model 
turns out to be a special case of the model considered 
here. The dynamical equation is particularly suited to 
capturing some essential features such as the diffusive 
spreading of the rain field through a fluid atmosphere. 
The underlying stochastic differential equation of the 
Yoo- Valdes-North (YVN) model includes a second time 
derivative term which may help in modeling the high 
frequency behavior of the rainfall spectrum. However, 
the model has more parameters and has so far been ap- 
plied only in the low wave number regime, and it is not 
yet clear how well the model describes rainfall statistics 
at the smaller spatial scales considered in this paper. 

The present model is an outgrowth of previous work 
by Bell [1987] and Bell et al [1990] who constructed 
a model to simulate the statistics of the gridded rain- 
fall data from GATE. They introduced a gaussian ran- 
dom field on a discrete 4 km grid whose Fourier am- 
plitudes were assumed to obey a stochastic dynamical 
equation similar to the one used here. A nonlinear map- 
ping was used to generate the gridded rain rate field 
with a rain-rate probability distribution matching the 
observed lognormal distribution. Unlike that work, the 
present model directly describes the stochastic evolu- 
tion of the instantaneous point rain rate field which can 
be averaged to any desired length and time scale. More- 
over, it yields the full space- time covariance function in 
terms of a small number of model parameters. Once 
they are determined from data, the model allows us to 
relate statistics at different scales, whereas the origi- 
nal model of Bell [1987] was applicable only to a single 
spatial grid scale and, moreover, the form of the spatial 
covariance function had to be specified as an explicit 
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function of spatial separation, in effect increasing the 
number of empirical parameters. However, there is a 
price to pay - the present model is restricted to dis- 
cussing only the second moment statistics. 

As we shall see in later sections, our model spectrum 
naturally leads to power law scaling behavior of the rain 
statistics at small spatial scales suggesting an underly- 
ing fractal structure of the rain field. Multifractal rain- 
fall models based on a somewhat different framework 
have been proposed by Schertzer and Lovejoy [1987] 
and Gupta and Waymire [1990] extending earlier ideas 
of Lovejoy [1982], Lovejoy and Mandelbrot [1985] and 
Waymire [1985]. These models describe the moments of 
a spatially averaged rain field in terms of multifractals, 
which treats the rain field fluctuations (i.e., deviation 
from the mean) as being governed by an underlying self- 
similar multiplicative random cascade process. 

The paper is organized as follows. Section 2 de- 
scribes the rainfall model. Here we collect a number 
of mathematical formulas related to the statistics of 
area-averaged rain rate field. Of particular interest to 
us is the way these statistics depend on the averag- 
ing length scale. In section 3 we describe the gridded 
TOGA COARE precipitation data set and the proce- 
dure for fitting the model to data and estimating the 
model parameters. Section 4 is devoted to a discussion 
of the results focusing, in particular, on the power law 
dependence of several rain statistics on the averaging 
length scale perhaps associated with fractal nature of 
the rain rate field. Some concluding remarks and direc- 
tions for future work are presented in the final section. 
Three appendices containing some mathematical details 
are included. 

2. Modeling Rainfall Variability 

In this section we describe the framework used in this 
paper to study rainfall variability. One of our main in- 
terests is to investigate how rain statistics depend on the 
averaging length and time scales. In this paper we con- 
centrate on the statistical properties of rain rate fields 
spatially averaged over various length scales. Gridded 
remote sensing precipitation data from ground based or 
satellite measurements are particularly convenient for 
such studies, since the gridding allows us to rescale the 
grid to any desired integral multiple of the original grid 
size at which the data is supplied, by simply grouping 
the grid boxes together into larger units. A number of 
basic statistics are calculated from the data for a range 
of grid sizes L and compared with model predictions. 
One of the main ingredients of the model is an ana- 


lytical form of the rainfall spectrum which is described 
next. 

2.1. A Spectral Model of Rain Rate Covariance 

In this subsection we give a brief account of the 
model. The underlying assumptions of the model as 
well as many of its salient features have already been 
described in detail in BK96. A set of results that are 
useful in fitting the model to data and checking some 
of its predictions are derived here. 

As was noted in Bell [1987] and Bell et al [1990], 
rain rate fields have the characteristic feature that the 
time scale over which spatial averages remain correlated 
depends on the size of the averaging area: the larger 
the averaging area, the longer the time over which the 
average rain rate remains correlated. This behavior, 
typical of turbulent fluids, must be captured by a re- 
alistic model of rainfall statistics. Although local rain 
statistics are expected to be affected by inhomogeneous 
and nonstationary directional effects, such as occur in 
rain bands, we are interested in a model description 
that is valid over a large enough region and sufficiently 
long time for such directional effects to average out, 
so that the rain statistics can be treated to a good ap- 
proximation as spatially homogeneous and isotropic and 
stationary in time. These assumptions imply that the 
space-time covariance function of the rain field at the 
points x, x' and times £,£' has the form 

c(x, t; x', t') = c(p, T ) = c(p, T ) 

where r = t* - t, p = x' — x and p = |p|. It can be ex- 
pressed as a Fourier integral 

c(p, r) = ^ 3/2 J dwj d 2 ke <(k P— U ' T) c(k,u>) , (5) 

where c(k,u;) represents the (unnormalized) rain rate 
power spectrum. Isotropy implies c(k,u;) = c(k,u) 
where k = |k|. The spatial Fourier amplitudes 

a(k,t) = J d 2 xe _lk ' x iJ / (x, t) 

of the field i?'(x,£) defined by Eq. (4) are assumed to 
obey the first order linear stochastic equation 1 

~a(k,t) = -^-a(k,t) + f{k,t) . (6) 

*A more general possibility is an equation of the form 

-^a(k,t) = - /d 2 k'— a(k',t) + /(k,t) , 
dt J T kk t 

which allows coupling among different Fourier modes. 
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Eq. (6) represents what is sometimes referred to as a 
first order autoregressive [AR(1)] process [Jenkins and 
Watts , 1968]. It is formally analogous to the familiar 
Langevin equation for the velocity of a particle execut- 
ing thermal Brownian motion in a viscous fluid. Here 
/(k,t) is a white noise force satisfying 

</(k, t)f* (k', 0) = (2 tt)F 0 5( k' - k) S(t' - 1) , (7) 

where in the right hand side we have introduced the 
Dirac <5-finetion, and the asterisk denotes complex con- 
jugation. In Eq. (6) r* represents the characteristic 
relaxation time of a(k, t) taken, as in BK96, to be of 
the form 

Tk = (1 + k 2 Ll) 1 + l/ ’ ^ 

where To and Lo are characteristic time and length 
scales. The length scale Lo effectively separates the 
long and the short wavelength regimes of the spatial 
fluctuations of the rain rate field. As a consequence of 
Eqs. (6) and (7) the power spectrum has the typical 
“red noise” form 


W>») = -2T-=5- ( 9 ) 

The analytic form of Tk in (8) is suggested by the form 
chosen by Bell [1987] and Bell et al [1990] who, in 
turn, were motivated by the empirical results of Laugh - 
lin [1981] from GATE data. As already noted, it con- 
tains the North-Nakamoto model as a special case with 
v — 0. It may be pointed out that in the original work 
of Bell [1987] and Bell et al. [1990], the form of the spa- 
tial covariance function was separately introduced from 
an empirical fit, whereas in the present model it is im- 
plicitly contained in the analytic form of t*.. Because of 
the fractional exponent, the stochastic dynamical equa- 
tion governing the evolution of the rain rate field now 
takes the form of an integro-differential equation rather 
than a pure differential equation as in the case of NN 
and YVN models. Physically, it can be thought of as 
representing a nonlocal “fractal” diffusion process by 
which rain moves through a turbulent atmosphere. In 
particular, the negative values of v , required to fit the 
data sets studied here appear to better describe the high 
spectral variability of rain rates at small spatial scales 
- an important factor in extrapolating statistics to rain 
gauge scales. 

The four model parameters F 0 , To, L 0 and v char- 
acterize the second moment statistics. No specific as- 
sumption about the actual probability distribution for 
the random variable i?(x, t) is introduced at this stage 


of the model applications. It should be noted how- 
ever that the model remains incomplete without such 
an assumption since the higher order moments cannot 
be represented within the model. Observationally, area 
averaged rain rates obey an approximately lognormal 
or gamma distribution, at least on smaller space and 
time scales. 

Substituting the explicit form of the model spectrum 
in (5) and performing the integral over w and the an- 
gular integral in the k-plane yields 

TOO 

c(p, t) = / dk kJo(kp)c(k,r) , (10) 

Jo 

where Jo{x) is the Bessel function of the first kind of 
order zero, and 


c(h,r) = y/n/2FoTke |r|/Tfc (11) 

represents the lagged covariance of a(k,t) defined by 
the equation 

(a(k, t)a*( k', t')) — (27 r)c(fc, r) 5(k' — k) , (12) 

and is the inverse temporal Fourier transform of c{k,u). 
At zero lag (t = 0) carrying out the remaining k- 
integration in (10), the spatial covariance of the point 
rain rate can be expressed in the simple closed form 

c(p, r = 0) = 7o C v (p/L 0 ) , (13) 


where, for convenience, 70 is defined by the relation 

Fo = y/2ftr(l + u)(Ll/Toho , (14) 

r(z) = f£°t z ~ 1 e~ t dt is the Euler P-function. and 


C v (z) = {z/2YK v {z) (15) 


is a function related to the modified Bessel function 
K v (z). For z <SC 1, C„(z) has the behavior 


1 1 / y \ 2l/ 

CM = ±I» + ^T(-u) (|) + 0(z 2+2v ) + 0(z 2 ) . 

(16) 

When v > 0, the first term dominates and C u {z) -+ ^r(^) 
as z — ► 0. In this case the variance of the point rain 
rate c(0, 0) = ^ / yoT(u) is finite. However, as we shall 
see, like the GATE data, the TOGA-COARE data re- 
quires v < 0 and consequently the covariance function 
exhibits a power law singularity 


c(p,0)~i7or(|H)(^/2io)- 2M (17) 
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as the separation p — ► 0. The origin of this singular 
behavior can be traced to the fc-dependence of the char- 
acteristic time scale r k governing the decay rate of the 
Fourier mode a(k, t ) ; for the long wavelength modes 
(k % 0) Tk approaches a constant value t 0 , whereas for 
the short wavelength modes (fc » 1 /Lo) Tfc behaves like 
£-2(i-M). Note that for values of the exponent v in 
the range —1 < v < 0 (as seems to be the case for 
both GATE and TOGA-COARE data sets), Tk 0 as 
k — > oo. For the NN model ( v = 0) the singularity in 
c(p, 0) weakens to become logarithmic. 

Next, we turn to the space-time covariance function 
of area-averaged rain rate, 

Caa'(s,t) = (ft A (t)K A ,(t + t)) 

where A, A! are two L x L boxes whose centers are 
separated by s and the prime on the rain rate variables 
denotes deviation from the mean, as in Eq. (4). Starting 
from the formula 

Caa' (s, t) = -^J^d 2 x j d 2 x' c(s + x' - x, r) (18) 

our spectral model yields at zero lag (r = 0) the expres- 
sion 

CaA'( s, 0 ) = 470 f f 1 d £ 2 (1 - fi)(l - £2) 

Jo Jo 

* Mk+ils)- !19) 

See BK96 Appendix A for details of the derivation. 
Upon setting s = 0 in Eq. (18), the variance of the 
area-averaged rain rate in a grid box A of area A = L 2 
is given by 

a\ = <7 2 (L) = Caa'{s = 0,r = 0) 

= 4ioG(v;L/Lo) , (20) 

where G(u ; z) is the double integral 

G(u; z) = f l fd^ (1 - ^)(1 - 6) 

Jo jo 

x (z\Jei+d) . ( 21 ) 

The power law singularity of the spatial covariance 

of the point rain rate as a function of separation in turn 
leads to a simple scaling behavior of the statistics of 
area averaged rain with the size of the area. Using the 
first two terms of the expansion (16) in Eq. (21), which 
dominate as long as — 1 < v < 0, we have 

G(y- z) « lr(-M) + \r(\v\)H(v)z~ 2 M , (22) 


where H(v) denotes the integral 


H(u)= f / 1 deid6(l-{i)(l-6)({? + £)•'. (23) 

JO Jo 

Hence, for small L (more precisely L < Lq), a 2 (L) can 
be expressed in the form 


a 2 (L ) «oo + boL 2 , 

(24) 

where ao and 6o are constants given by 


ao = ^7dT(-M) , 

(25) 

bo = 27 oT(\v\)H(v)(2L 0 ) 2M . 

(26) 


Thus the model predicts that as L — ► 0, the variance 
of area-averaged rain rate cr 2 (L) exhibits a power law 
singularity: 

a 2 (L) - . 

We note that although the point variance is infinite in 
the model, it is not necessarily inconsistent with experi- 
ence since observations always involve smoothing of the 
point rain rate field in space and/or time, and variances 
of smoothed rain field are finite in the model. 

The lagged correlation function 

®AA' (s ; r) = ('AA'is, t) (27) 

a A 

is a rapidly decreasing function of the separation s. The 
model correlation function is markedly nonexponential. 
Its fall-off rate is conveniently characterized by the “in- 
tegral correlation length” Ai n t (L) [which would be the 
(l/e)-folding distance for an exponentially decaying cor- 
relation function] defined in the isotropic case through 
the relation 

Af„t (L) = [°°ds s$ AA .( s, r = 0) . (28) 

Jo 

It is related to the model spectrum through 

A ‘ nt(L) = ^Z) c(k = 0 ’ r = 0) ’ (29) 

which upon using Eq.(20) reduces to 

AintW = \t(1 + v)Ll/G{v, L/Lo) . (30) 

From Eq. (22) it then follows that Aj nt (T) vanishes as 
l) v \ as L -+ 0. 
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We can also define an integral autocorrelation time 
Tint^) (which would be the 1 /e-folding time for an ex- 
ponentially decreasing autocorrelation function) through 
the formula 


TOO 

Tint(L) = dr $aa(s = 0, r) . 
Jo 

It can be reexpressed in the form 


(31) 


(32) 


where C^uj) denotes the Fourier transform of the lagged 
autocovariance function Caa(& = 0, r). Evaluation of 
the integral in the numerator (see Appendix A) yields 


Ti n t(L) = Tq 


T(l + ^) G(1 4- 2i/; L/Lq) 
T(2 + 2u) G(u;L/Lo) 


(33) 


The values of v that best describe rain data in fact ap- 
pear to lie in the narrower range — 1 /2 < < 0 (rather 

than the range — 1 < v < 0 considered above). The nu- 
merator of Eq. (33) is then finite and consequently the 
model predicts that T int (L) — > 0 as as L — ► 0. In 
other words, the “time scale” for fluctuations of average 
rain rate in a box area L 2 gets smaller with decreasing 
box size. 


The case v — — 1 presents a particularly simple 
and analytically tractable special case of our spectral 
model. In this case, the characteristic length Lo drops 
out of the model. All the spatial Fourier modes ex- 
cited by the random white noise source decay expo- 
nentially with the same characteristic e-folding time To- 
The Fourier transform of the spatial covariance c(p,0), 
namely, c(k, 0) is a constant, which in turn implies that 
the point rain rate field at equal times is ^-correlated: 

c(p,0) oc S(p). 


2.2. Some Other Rain Statistics 

The spectral model presented above makes a number 
of predictions about the change of rain statistics with 
scale. With a few additional assumptions these predic- 
tions can be shown to have interesting consequences for 
the asymptotic behavior of the “conditional” statistics 
of rain as well. We now turn to a brief description of 
these statistics. 

Consider remote sensing measurements of precipita- 
tion over a given area and time period and the resulting 
rain maps gridded on an L x L square grid. The ensem- 
ble mean rain rate (R) is estimated as the average of all 
the data within the data set, and, due to the assump- 
tion of statistical homogeneity, is independent of the 
resolution scale L. However, other statistical properties 
of the field, such as the fraction of grid boxes contain- 
ing nonzero rain (“the rainy area fraction”) p(L) and 
the variance of the grid box-averaged rain rate <t 2 (L) 
are expected to depend nontrivially on L. Two other 
related statistics are the mean r c (L) and the variance 
o 2 (L) of the grid box- averaged rain rate conditional on 
the box containing non zero rain. It is easy to show that 
(see, e.g., BKK, Appendix A) these are related through 


(R)=p(L)r c (L), 

v 2 {L) = p(L)[<rl(L) + r c (L) 2 } - ( R ) 2 
It is convenient to define the dimensionless ratio 

M = p{L) = a c (L)/r c (L). 

Using Eq. (35), Eq. (36) can then be written as 
cr 2 (L) 1 + ( 1 ? 


1 + 


(R)* 


P{L) ' 


(35) 

(36) 

(37) 


(38) 


In this case the covariance function of the area-averaged 
rain rate at zero lag can be easily computed for v = -1 
to yield the explicit formula 

- (^X*- ¥)(*-¥) 

x w (x) w (^) ’ (34) 

where W (a;) is the window function 

Wix) — / 1 < 1 

o \ x \ > 1 * 

The correlation between the rain rates in disjoint grid 
boxes vanish: the integral correlation length Af nt (L) is 
simply equal to the box size L. The variance of area- 
averaged rain rate a 2 (L) is proportional to 1/L 2 . 


Based on analysis of rain gauge data from Darwin, 
Australia and Melbourne, Florida and radar data from 
GATE, Short et al [1993] have suggested that p is rel- 
atively constant over a range of length scales, averag- 
ing times, types of data (radar or rain gauge) and rain 
climates. In particular, p is insensitive to the averag- 
ing length scale L. This was also corroborated by the 
results of BKK from precipitation data obtained from 
SSM/I instruments on board NO A A DMSP satellites 
F10 and Fll. 

If the quantity p(L) remains finite as L — ► 0, then 
Eq. (38) together with the asymptotic behavior (24) 
predicts that p(L) would vanish as L 2|l/| in this limit. 
The limiting behavior p(L) — > 0 as L — > 0 in turn im- 
plies that the local instantaneous rain rate field i?(x, t) 
vanishes almost everywhere , i.e., everywhere except a 
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set of measure zero - a familiar property of spatial frac- 
tals. We shall return to this in greater detail in section 
4 when we present the results of our data analysis. 

3. Testing the Model Predictions with 
Ground based Radar Data 

3.1. The TOGA COARE dataset 

The radar data we use to test the spectral model 
described above were collected during TOGA COARE 
[Webster and Lukas , 1992] carried out in the West- 
ern Pacific Equatorial Warm Pool during the period 
November 1992 to February 1993. The measurements 
were made with a pair of calibrated Doppler radars (la- 
beled TOGA and MIT) on board two research ships lo- 
cated near the center of the Intensive Flux Array (IFA) 
- a region in the shape of a quadrilateral, roughly cen- 
tered at 2° S, 156° E. Radar returns weaken with dis- 
tance from the radar, and quantitatively reliable rain 
rate estimates were provided out to a distance 145 km 
from each radar. The fields of view (FOVs) of the two 
radars overlapped substantially when they were both 
operating. The measured radar reflectances were con- 
verted into average surface rain rates on a 2 km spatial 
grid using an empirical Z — R relation. The gridded 
rainfall data set used in the present work is described 
by Short et al. [1997]. The total 101-day period of ob- 
servation was divided into three “monthly” cruises cov- 
ering the periods November 11 to December 10 1992, 
December 15 1992 to January 18 1993 and January 23 
to February 23 1993 designated respectively as Cruises 
1, 2 and 3. Rain maps were available roughly at in- 
tervals of ten minutes with occasional long gaps during 
which data from at least one of the radars were missing. 

3.2. Statistical Analysis of the Data 

The gridded data set obtained from the TOGA 
COARE observations yields estimates of the various 
area-averaged rain statistics introduced previously, al- 
lowing one to test the model predictions. We now briefly 
describe how these statistics are obtained from the data. 

The statistics are derived using rain estimates from 
square areas 128 x 128 km 2 centered at the location 
of each of the two ships. The locations of the areas 
changed slowly with time because of the slight drift of 
the ships during the month-long experiments. However, 
the effect of this on the statistics, if any, is expected to 
be negligible. Each area is subdivided into boxes Lx L 
in size, with L — 2 n km (n = 1,2, . . . ,7), and the 
quantities ( R ), <j 2 (L), p{L ), r c (L) and cr c (L) defined 


in the previous section are computed for each L from 
the data. Only those Lx L boxes in a rain map were 
used for which at least 95% of the box had valid data. 
This was necessary since the radar images frequently 
had data missing along a radial line emanating from the 
center of the circular FOV where the radar is located. 
This data dropout phenomenon affected the grid boxes 
located close to the origin through which this radial 
line passed, especially at the smaller grid length scales 
L = 2, 4, 8 and 16 km. 

For comparison purposes, we also computed the same 
spatial statistics for the GATE Phase I (1716 images) 
and Phase II (1512 images) data from the 280 km box 
centered at the location of the ship. For GATE data 
the basic grid size is 4 km. 

We determined the integral correlation time T[ nt {L) 
for each L as follows. Mean rain rate time series were 
constructed for each Lx L box within the selected ar- 
eas for the three cruises. Again only boxes with at least 
95% of the area containing data were included in the 
time series in order to contend with the data loss phe- 
nomenon mentioned above. The lagged autocovariance 
function Caa{®, r) averaged over all boxes in the chosen 
areas was computed for each L at lags r of multiples of 
10 min. 

As an example, the lagged autocorrelation function 
$l(t) = $A4(0,r) 

for L = 128 km for the six data sets are shown in 
Figure 1. For each data set a numerical estimate of the 
time integral r mt (L) = J 0 °° 4> L (r)dr was obtained. This 
is discussed in more detail later. 

The following subsection describes how the parame- 
ters 7 q, L 0 , v and r 0 are determined by fitting the data 
to the statistics of gridded averages predicted by the 
model. 

3.3. Estimation of the Model Parameters 

The parameters 70, io and v determine the spatial 
statistics at zero lag. They are estimated by fitting the 
variance of the gridded rain rate average a 2 (L) to the 
asymptotic formula (24) valid for small L. The fit deter- 
mines the exponent v and the coefficients a 0 and 60 • Eq. 
(25) then determines the parameter 70 and Eq. (26) fi- 
nally determines the characteristic length scale Lo- The 
procedure was carried out for the TOGA COARE as 
well as both GATE Phase I and Phase II data sets. 

An estimate of the characteristic time To for each 
data set can be obtained by averaging over the values 
individually computed by solving (33) for r 0 using the 
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Lagged Autocorrelation for TOGA CO ARE Data Sets 



Figure 1 . Plot of lagged autocorrelation functions 
for L = 128 km square boxes centered at the 
location of the TOGA and MIT radars for each of the 
three cruises. 

values of the integral correlation time r\ nt (L) for various 
box sizes L and the previously determined parameter 
values of v and Lq . Since there is only a finite length T 
of data available, one cannot actually calculate T| n t(£) 
from the data itself. The data can only supply an esti- 
mate of (31) with a finite cut-off, 

/*Tmax 

T int(L, TO] Tmax) — / $ L {r)dT (39) 

Jo 

with r max less than T. In principle, we could solve (39) 
for to using the model from the data-derived value of 
fi n t for a chosen cut-off r max . In practice this is dif- 
ficult, since the model-predicted r 0 -dependence of the 
function T int (L, To; r max ) is nonlinear and complex to 
calculate. Instead, we adopt an iterative method which 
we describe next. 

There is an uncertainty inherent in choosing T max in 
order to get the “best estimate” of tq. It arises because 
the tail of the autocorrelation function suffers from a 
combined effect of statistical noise and possible secular 
variations of the rainfall pattern over longer time scales 
not represented in our spectral model. The choice of the 
cut-off is guided by several considerations. On the one 
hand, r max should be large enough so that the “resid- 
ual” effect of the tail in the integral (31) is small. Some 
guidance on an acceptable value of r max is provided by 
the model’s prediction of the the truncated integral. 
However, the empirical lagged autocorrelation suffers 


from sampling error because it is obtained from a finite 
time series of length T « 30 days. The cut-off r max 
must therefore satisfy the condition r max <C T. Also 
r max should be chosen to be small enough so that the 
effect of slow secular modulations of the lagged autocor- 
relation function caused by physical phenomena such as 
a diurnal cycle or various atmospheric waves, which are 
not accounted for in our simple model spectrum, is min- 
imized. Evidence of such longer time scales was most 
obvious in the Cruise 2 MIT and the Cruise 3 TOGA 
radar data sets. These modulations can probably be at- 
tributed to various large scale atmospheric waves known 
to influence rainfall in the equatorial Pacific. We shall 
revisit some of these issues in the next section. Keep- 
ing these factors in mind, the cut-off parameter r max 
employed in our computations was set at 12 hrs. 

The choice of the upper limit of integration is some- 
what subjective. The tails of the observed correlation 
functions were substantial for several of the data sets, 
most notably, TOGA Cruise 3, for which it extended 
well beyond 20 hours, at which lags one might expect 
effects of diurnal variation and atmospheric waves to 
influence the statistics. We therefore chose a relatively 
small cut-off r max = 12 hours with the assumption that 
the spectral model accurately describes the observed 
correlations up to the lag r = r max . 

Next, we turn to the problem of estimation of r 0 . 
Using (39), Eq. (31) for the integral correlation time 
r in t(L) = 7int(^ 5 To) can be expressed as the sum 

roo 

Tint{LyTo) = Ti n t(T,To; Tmax ) + / $L(r)dr (40) 

J Tnax 

If the left-hand side of (40) were known, Eq. (33) would 
immediately yield r 0 . However, while the first term 
on the right-hand side of (40) is obtained from data, 
the second term can be calculated from the model only 
when to is known. We therefore solve Eq. (40) iter- 
atively: to begin the process, we obtain a first guess 
for To by ignoring the contribution of the integral; for 
all the following steps of the iteration we evaluate the 
integral on the right-hand side using the model $i(t) 
with the value of To obtained from the previous step. 
For ease of computation, we replace the square box by 
a circular disk of equal area. This approximation is suf- 
ficient for our purposes. Appendix B gives some details 
of the model calculations for a circular area. In prac- 
tice, for the TOGA COARE data, the correction to to 
due to the introduction of a cutoff is of order 5%, and 
it was unnecessary to go beyond the first iteration. 

The iterative method applied to GATE Phase I data 
with T max chosen to be 18 h gives an estimate of tq 
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within about 1% of what was reported in BK96 from a 
direct fit to the lagged autocorrelation function. 

4. Results and Discussion 

The full set of parameter values for which the model 
best fits the TOGA COARE data are given in Table 1, 
along with the corresponding values for GATE Phase I 
determined in BK96 using a somewhat less exact (and 
more arduous) trial-and-error curve fitting method in- 
stead of the analytical method adopted here. 

4.1. Model Parameters v , Lq and 70 

Table 1 summarizes the TOGA COARE model pa- 
rameters. We see that the values of the power law expo- 
nent v range between about —0.21 to -0.34 indicating 
a much stronger singular behavior than was found for 
the GATE Phase I data ( 1 / — -0.11). The character- 
istic time scales To lie in the range of roughly 4.5 to 
8.2 hours, considerably shorter than the GATE Phase 
I value To « 13 hours. The characteristic length Lq 
ranges between about 54 and 94 km indicating a some- 
what sharper fall-off of the spatial correlations com- 
pared to GATE Phase I for which Lo — 104 km. The 
results exhibiting the quality of the fit are shown in Fig- 
ure 2 for the TOGA and MIT radars. It is seen that in 
most cases the fit is quite good up to grid size L — 64 
km, beyond which the approximation on which Eq. (24) 
is based is expected to break down. Our results seem to 
be consistent with the model prediction that the vari- 
ance cr 2 (L) has a power law behavior cr 2 (L) ~ Z,~ 2 M as 
L — > 0. Robustness of the parameterization was checked 
by comparing the variance of area-averaged rain rate 
(<7 2 ) t h for an area A = L 2 with L = 128 km predicted 
by the exact formula (20), rather than the asymptotic 
formula (24), against the observed variance (a^) 0 bs for 
the 128 km box centered at the two ship locations. Ta- 
ble 2 shows the comparisons, which Eire quite good. 

For GATE data we find that for L between 4 km and 
56 km, the variance cr 2 (L) is quite accurately fitted (not 
shown here) by the formulas 

2/ t\ f 17.3L' 0 - 24 - 4.66 mm 2 h~ 2 Phase I 

K ’ \ 15.0L- 0 38 - 1.92 mm 2 h -2 Phase II ’ 

which through Eqs. (24) to (26) yield the parameter 
values 

7o = 1.02 mm 2 h" 2 , v = -0.12, L 0 = 93.8 km 
for Phase I and 

70 = 0.63 mm 2 h“ 2 , v = -0.19, L 0 = 82.7 km 


Variance as function of area for TOGA COARE 




Figure 2. Comparison between the observed variance 
of area-averaged rain rate a 2 (L) for the TOGA COARE 
data and the corresponding spectral model predictions 
for various box sizes L ranging from 2 km to 128 km 
using Eqs. (24)-(26) and the parameter values listed in 
Table 1. 
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Table 1. Parameter Values for Rain-Rate Covariance Model. 


Dataset 

7o (mm 2 h~ 2 ) 

V 

L 0 (km) 

T 0 (h) 

( R ) (mmh -1 ) 

GATE Phase I 

1.0 

-0.11 

104. 

13.0 

0.492 

TOGA Cruise 1 

0.067 

-0.335 

94.06 

6.8 

0.139 

MIT Cruise 1 

0.086 

-0.297 

73.89 

5.8 

0.134 

TOGA Cruise 2 

0.616 

-0.239 

53.81 

5.0 

0.351 

MIT Cruise 2 

0.206 

-0.205 

70.40 

8.2 

0.229 

TOGA Cruise 3 

0.127 

-0.290 

61.04 

5.2 

0.155 

MIT Cruise 3 

0.180 

-0.259 

64.94 

4.5 

0.200 


Parameters for the model spectrum obtained from fits to radar data from Phase I of GATE and from the 
ships MIT and TOGA during the three TOGA COARE cruises. Average rain rate for each dataset is given in 
the last column. 


for Phase II. [The original GATE Phase I parameters 
listed in Table 1 lead to the asymptotic formula cr 2 {L) = 
16.80L- 0 * 22 - 4.89 mm 2 h“ 2 (see BK96 Appendix A), 
whose coefficients and the power law exponent are well 
within the 95% confidence interval of the corresponding 
quantities in the empirical least squares fit given above.] 

We also explored the dependence of the model pa- 
rameters on the mean rain rate (ii). While no system- 
atic dependence could be discerned for L 0 , tq and v, 
the strength parameter 7 0 was found to depend on ( R ) 
according to the simple formula 

70 = k(R) 2 , (41) 

where k is a dimensionless constant. The fit in Figure 
3, which also includes points corresponding to GATE 
Phase I «J?> = 0.492 mmh- 1 ) and GATE Phase II 
(( R ) — 0.386 rrirnh" 1 ), yields a value k » 4.33. 

4.2. Model Parameter tq 

Next, we turn to a discussion of the dependence of 
the characteristic time scale r\ n t(L) on the spatial av- 
eraging scale L. The plots in Fig. 4 show reason- 
able agreement between the values of r[ nt (L) determined 
from data and the corresponding theoretical predictions 
computed for various box sizes L using the model pa- 
rameters. Linearity of the log-log plot for small L con- 
firms the power law dependence Ti nt (L) ~ ( L/Lq ) 2fzy| 
as L — ► 0 expected from our model. In BK96 it was 
found that the integral correlation time estimates for 
the GATE Phase I is described quite well by the power 
law form r int (L) = 0.67L 0 * 49 h. Again the exponent of 
the power law is rather different from the expected value 
0.22 obtained from the spectral model fit. However, it 


Plot of y 0 vs. Mean Rain Rate 



Figure 3. Plot of the parameter 70 against the mean 
rain rate ( R ) for the TOGA COARE and GATE data 
sets. The continuous curve shows the least squares fit 
to the formula (43) 
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Table 2, Comparison between theoretical and observed values of u\ for L = 128 km. 


Dataset 

(oi)th (mm 2 h 2 ) 

(^)obs (mm 2 h 2 ) 

TOGA Cruise 1 

0.107 

0.116 

MIT Cruise 1 

0.093 

0.105 

TOGA Cruise 2 

0.399 

0.371 

MIT Cruise 2 

0.176 

0.174 

TOGA Cruise 3 

0.107 

0.119 

MIT Cruise 3 

0.155 

0.169 


Comparison between values of the variance of area-averaged rain rate <j\ computed from the model using 
Eq. (20) and the parameters listed in Table 1 for L = 128 km and from the data for two 128 km boxes centered 
at each of the two ship locations. 


was also found there that while the model captures the 
overall decay of the autocorrelation function $l(t), it 
substantially underestimates the correlation at smaller 
lags, indicating a departure from the model spectrum 
for a first order process. Indeed, it was shown by Bell 
and Reid [1993] that this behavior could be accounted 
for by going over to a second order model, which of 
course introduces more adjustable parameters. 

A number of caveats in the estimates of ro should be 
pointed out. A fundamental uncertainty in the value of 
r int (L) arises from the fact that it must be estimated 
from a time series of finite length T which in our case 
is of the order 30 days. It is known [Bell, 1980] that 
if the true autocorrelation function of a normally dis- 
tributed random variable has an exponentially decaying 
form = exp(— r/r int ) then the integral correla- 

tion Tint? defined by Eq. (38) with a finite cut-off r m8L x, 
has variance 

<7 ['Hnt] = 4(r max /T)fi^ t , (42) 

where it is assumed that ri nt (X) < r max <T. A proof 
of Eq. (42) is given in Appendix C. Thus, in order to 
reduce the statistical uncertainty in the estimate, r max 
should be chosen to be as small as possible but large 
enough to include the entire range over which 4 >l(t) 
differs appreciably from zero. In our case, if we choose 
T ma x « 12 hrs, then it follows that the 2<r uncertainty in 
the estimate of r mt (L) obtained from a single time series 
could, in fact, be as large as 0.5fj nt . It is to be noted 
that, strictly speaking, this estimate holds only for a 
gaussian random variable undergoing a linear stochas- 
tic dynamical process characterized by an exponentially 
decaying autocorrelation function. However, the proba- 
bility distribution of rain rate variable is markedly non- 


gaussian and the autocorrelation functions possess a 
much more slowly decaying tail. These factors can lead 
to an uncertainty larger than what is predicted by Eq. 
(42). On the other hand, for the smaller grid boxes of 
size L = 2 n km with n < 6, the actual uncertainty 
should be considerably lower than the estimate (42) 
since the autocovariances were obtained by averaging 
over the results for the time series for the 4 7_n boxes. 
The time series are, however, not completely indepen- 
dent of one another, because of the spatial correlation 
of the rain rate in a box with others in its vicinity, and 
an exact estimate of the uncertainty would be difficult 
to obtain. 

4.3. Conditional Statistics and Spectral Model 

The spectral model described in section 2 allows one 
to quantify how the statistics of area-averaged rain rate 
depends on the averaging scale L. These statistical 
quantities, such as the variance <j 2 (L) and the integral 
correlation length and time scales, pertain to variabil- 
ity of the rain field including its zeroes, i.e., the re- 
gion where it is not actually raining. However, unlike 
most random variables in nature, the instantaneous rain 
rate has a finite nonzero probability of having the sharp 
value zero, in accordance with what is a matter of com- 
mon experience - namely that it is not raining in most 
places most of the time. Mathematically this is ex- 
pressed by the property that the probability distribu- 
tion function for the area-averaged instantaneous rain 
rate 7^4 (£) defined by (1) has a 5-function peak at zero 
followed by an (approximately lognormal) continuum 
[Kedem et al , 1990]. The probability p(L) for the value 
7Z A ^ 0 is expected to have nontrivial L-dependence 
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Plot of integral correlation time vs. L for TOGA COARE 




Box size L(km) 


Figure 4. Plot of the integral correlation time fi nt 
evaluated with a finite cut-off r max for the three TOGA 
and MIT cruises and various box sizes L ranging from 
2 km to 128 km. The curves represent interpolation of 
the values computed from the spectral model. 


Test of constancy of \i 



Figure 5, Scatterplot of 1 + a 2 (L)/(R) 2 vs. 1 /p{L) 
for all the six TOGA COARE data sets for averaging 
length scales L ranging from 2 km to 128 km (6 x 7 = 42 
data points). Least squares fit to a straight line through 
the origin yields the slope 1 + /x 2 « 10.2, i.e. p « 3.03 
for the entire TOGA COARE data set. 


and is an important statistic of the gridded rain field. 
As discussed in section 2, an additional property of the 
conditional statistics of rain, namely the constancy of 
the dimensionless ratio /i, allows one to draw some infer- 
ences about the asymptotic L-dependence of the func- 
tion p(L). We now turn to a discussion of this statistic 
for the TOGA COARE data set. It is, however, worth 
pointing out that in doing so we are exploring terri- 
tory that lies outside the scope of the original spectral 
model. 

Our data set provides a simple way to test the con- 
stancy of the ratio p [defined by Eq. (37)] as suggested 
by Short et al [1993]. Eq. (38) implies that the plot 
of 1 + o 2 (L)/(R) 2 vs. l/p(L) would be a straight line 
through the origin if p is a constant independent of L. 
A scatterplot of all the data points from the six data 
sets (2 ships, 3 cruises) and each value of L is shown 
in Figure 5. The linear fit seems satisfactory; the slope 
yields a mean value of p & 3.03 collectively for all the 
six data sets. However, there is a pronounced system- 
atic deviation from the linear fit in the “large* p(L) 
regime [i.e., p(L) ~ 1], which corresponds to large L. 
One can also evaluate p separately for each data set. 
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The individual values were found to be within a narrow 
range, with a low value of 2.63 for MIT Cruise 2 and 
a high value of 3.39 for TOGA Cruise 2. As noted in 
section 2, constancy of p (or more generally, finiteness 
of p as L — > 0) together with (24) leads to another im- 
portant conclusion - a power law scaling behavior of 
the rain probability p(L) at small L. Combining the 
formulas (24)-(26) and (41) with Eq. (38) we infer that 
for L Lq , 


l/p(I)ocl + /3(L/2Lo)- 2|,/1 , (43) 


where 


2kT(\v\)H{v) 

1 + K(-M)’ 


(44) 


H(v) is given by Eq. (23) and k is the dimensionless 
number introduced in (41). This immediately implies 
the power law behavior 



as L — > 0, which can be independently checked from the 
data. Fits to a power law p(L) oc L v for each of the six 
data sets are shown in Figure 6. A straight line fit on 
a plot of logp(L) vs. log L provides evidence of a power 
law behavior between L = 2 km and 64 km. In Table 
3 the observed exponent r/ obs is compared against the 
theoretically expected value r? t h = 2\v\ for each data 
set. We see that, although there is generally speaking a 
close agreement between the two sets of values, for four 
of the six data sets r? 0 bs is somewhat larger than 7/ t h 
contrary to what one would infer from Eq. (38). How- 
ever, it should be noted that according to the model 
the power law dependence on L becomes accurate only 
in the scaling regime L Lq. Consequently, it is more 
appropriate to compare r? tb with the power law expo- 
nent determined from p(L) in the small L limit. Indeed, 
a closer inspection of the graphs of p(L) shows a slight 
but systematic deviation at small L from the power law 
obtained by least squares fit, indicating a smaller slope 
there. In Table 3 we have also included an estimated 
exponent rj^ bs evaluated from the slope of a linear in- 
terpolation between L — 2 km and 8 km, which clearly 
improves agreement between the model and observa- 
tion. One can also see in retrospect that the constancy 
of p which gave rise to the power law scaling behavior 
of p(L ), itself becomes progressively more accurate in 
the limit of small L . 

For the GATE data we find that within the range of 
L between 4 km and 56 km the rain probability p(L) 
(L in km) is represented quite well by the power law 


Plot of p(L) vs. L for TOGA COARE 




Figure 6. Power law fit to the fractional rainy area 
p(L) for the three TOGA and MIT cruises and various 
box sizes L ranging from 2 km to 64 km. The exponents 
are listed in Table 3. 
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Table 3. Comparison between theoretical and observed values of the exponent 77 . 


Dataset 

Tjobs 

^7obs 

Vth - 2 M 

TOGA Cruise 1 

0.760 ± 0.026 

0.571 

0.669 ± 0.025 

MIT Cruise 1 

0.728 ± 0.025 

0.618 

0.595 ± 0.047 

TOGA Cruise 2 

0.464 ± 0.011 

0.382 

0.477 ± 0.057 

MIT Cruise 2 

0.540 ± 0.014 

0.417 

0.409 ± 0.054 

TOGA Cruise 3 

0.582 ± 0.015 

0.448 

0.581 ± 0.056 

MIT Cruise 3 

0.573 ± 0.015 

0.448 

0.519 ± 0.063 


Comparison between value of the power law index rj defined by the stipulated behavior p(L) oc L v of the 
fractional rainy area p(L) as function of the averaging length scale L determined from the data and the value 
predicted by the model. The second column gives the “mean” exponent obtained by an overall least squares fit, 
while the third column gives the exponent determined from the slope at small L. In the last column specifying 
the “theoretical” value 77 th = 2|i/| determined from the a 2 (L) vs. L plots [Figure 2 via the fit to (24)], we also 
include the error bars not reported in Table 1 . 


formulas 

m f 0.122(L/4 ) 0 * 61 Phase I 
P[ } \ 0.098(L/4)°* 68 Phase II ’ 

which are consistent with the expected behavior p(L) — > 

0 as L — y 0. [See also Ha et al. , 2002, who use a linear 
regression fit to describe p(L) as function of L], The 
exponents are, however, in each case quite different from 
the values suggested by the spectral model fit. We also 
find that for GATE, the quantity p(L) actually depends 
significantly on L. In fact, p(L) changes from about 
1.77 at L = 4 km to about 2.16 at L = 56 km for Phase 

1 and from about 1.89 at L = 4 km to about 2.22 at 
L — 56 km for Phase II. 

Like in the TOGA COARE case, the goodness of a 
simple power law fit is probably somewhat fortuitous 
for GATE as well, since deviations from a power law 
behavior is to be expected on the basis of the model, be- 
cause of the presence of “subdominant” terms in 1 jp{L) 
such as the constant term in (43) or other possible In- 
dependent terms arising from a weak L-dependence of 
p(L). Upon attempting to fit the data for p(L) to a 
form like (43) we find that the best fit values of the 
exponent 77 are slightly smaller than the values given 
above. However, they still differ substantially from the 
model-inspired value 2\u\. We have been unable to fully 
reconcile the apparent discrepancy between the two sets 
of exponents along the lines discussed above in connec- 
tion with the TOGA COARE data. 

One particular aspect that we did not address in this 
paper is the dependence of p(L) and other rain statistics 


on the value of the rain threshold, i.e., the smallest value 
of rain rate that is experimentally measurable. This is 
an important issue, since different instruments have dif- 
ferent thresholds and hence threshold dependence can 
greatly complicate intercomparison efforts for rain rate 
averages, especially conditional averages. 

4.4. Discussion 

From Table 1 it is evident that the model param- 
eters vary considerably among the data sets. These 
variations are almost certainly due in part to real sec- 
ular variations in the rainfall pattern. As discussed 
by Short et al [1997], the TOGA COARE region of 
the tropical Western Pacific experiences planetary scale 
equatorial waves including an annual cycle, intrasea- 
sonal oscillations such as Madden- Julian oscillations 
[Lau and Chan , 1985; Madden , 1986; Nakazawa , 1995] 
and westerly wind bursts related to El Nino-Southern 
Oscillation (ENSO) events. The first half of Cruise 2 
had an intensely active period associated with intrasea- 
sonal oscillation accompanied by passage of cloud su- 
perclusters and heavy rainfall. The rainfall time se- 
ries during this period also shows a quasi 2-day oscil- 
lation of the type described by Takayabu [1994] and 
Takayabu et al [1996]. They attribute this oscillation 
to westward-propagating inertio-gravity waves. By con- 
trast, the spatial rainfall pattern during much of Cruise 
1 appears to be much more homogeneous and isotropic, 
characteristic of small-scale convective systems. For a 
more detailed discussion see Short et al [1997]. 
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5. Conclusion 


In this paper we have presented a model of secopd 
moment statistics of rainfall based on a linear stochas- 
tic dynamical equation. The model makes definite pre- 
dictions about the statistical behavior of area-averaged 
rain rate in both space and time. In particular, the 
model predicts simple power law dependence of the rain 
statistics on the averaging area as the area tends to 
zero, perhaps reflecting a fractal nature of the rain field 
at smaller spatial scales. These predictions are tested 
with a radar-derived gridded tropical oceanic precipita- 
tion data set and the model is shown to describe the 
data quite well within the intrinsic limitations of the 
model assumption of statistical homogeneity, isotropy 
and stationarity. In particular, the intimate relation- 
ship between the dependence of the correlation length 
and time scales of the rain field on the spatial smooth- 
ing predicted by the model is in accordance with data. 
It should, however, be emphasized that in extrapolating 
the model-derived limiting behavior of the rain statis- 
tics deduced from a data set on a finite grid as the aver- 
aging scale goes to zero, we are effectively probing the 
“subgrid” structure of the rain field, i.e., the structure 
on scale finer than allowed by the instrumental resolu- 
tion. 

To conclude, we briefly mention some directions for 
future work. Like the spatially averaged rain rate, the 
spectral model also predicts relationships among the 
power law scaling exponents for the time-averaged rain 
rate as well. It would be interesting to see how well the 
model is capable of describing statistics of rain gauge 
data averaged over various time intervals. In a sequel to 
this paper, we will investigate the statistical properties 
of time-averaged local rain rates as measured by rain 
gauges using our model spectrum. Our model provides 
a convenient framework for addressing intercomparison 
problems that involve averaging the precipitation field 
over many different length and time scales. Recently the 
model has been used to investigate the problem of val- 
idation of satellite measurements of rainfall using rain 
gauge measurements on the ground [Bell and Kundu , 
2002], Finally, simple power law behavior of the various 
statistics of rain fields with respect to rescaling of the 
spatial averaging grid at small scales suggests a model- 
independent explanation based entirely on the scaling 
properties of the various statistical quantities involved. 
We plan to explore this possibility in future publica- 
tions. 


Appendix A: A Closed Analytic 
Formula for r int (L) 

In this appendix we give a derivation of the formula 
(33) for the integral correlation time Ti nt (L) and de- 
scribe its limiting behavior as L — > 0. 

Equation (20) for the variance of the area-averaged 
rain rate in a grid box A of area A = L 2 gives the 
expression 

°a = 47oG(i'; L/L 0 ) , 

where G(v;z) is defined as the double integral (21). 
On the other hand, the lagged autocovariance function 
C AA (0, T ) of rain rate averaged over the same area is 
given by (BK96, Appendix A) 

Caa(0,t) = \j™j™ dk idk 2 Q 2 (^Yj 

(^r) c ( k ’ T ) 

- vf *r £**•*(¥) 

x« 2 (^) 


where Q(x) = sin 2 z/:r 2 is the Bartlett filter function. 
Setting t = 0, we have an alternative integral represen- 
tation of : 




■ ?rr— «•(¥)«•(¥) 

l~2 /»00 r c 

V 7T Jo Jo 


Tk 


' dkidki G 2 (kiL/2) Q 2 (k 2 L/ 2) 

[1 + (k 2 + k 2 )L 2 }^ ’ 

(Al) 


where in the last step we have made use of the explicit 
form of Tfe given by (8). Comparing the two represen- 
tations of o 2 a we obtain the useful integral identity 


n °°dkidk 2 Q 2 (kiL/2) g 2 (JfejL/2) 

[1 + (k 2 + k 2 )L 2 ]^ 

- f^m G(v '’ LIL ° h <A2) 


Next, consider the formula for the Fourier transform 
of C AA (0,t) : 


C'a M = I J™ J~dkidk 2 g 2 (^) 
*g 2 (Jy) ■ 
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Using the explicit form of the model spectrum (9) with 
uj = 0, we get 

c A (u = 0 ) = | Fo J~J~dhdk 2 g 2 ) 

[2 r, _2 r rdhdk 2 g 2 ihL/2) g 2 ( k 2 L/2 ) 

Vtt or °y 0 y 0 [l+^+^ig]^) • 

In view of the identity (A2) this reduces to 

Ca(w = 0) = 4 ^/~ 7° T 0 p(2 + 2z/; ’ 

(A3) 

where we have also used (14) to substitute for Fo in 
terms of the standard model parameters 70, Lq and v. 
Using Eqs. (20), (A3) and (32) we obtain Eq. (33). 
The integrals were numerically evaluated with the help 
of Mathematica [v. 4; see Wolfram , 1999]. 

Finally, we examine the limiting behavior of T[ Rt (L) 
as L — ► 0. We note that for the range of values of 
the exponent v of interest to us, -1/2 < v < 0, (or 
0 < l+2iv < 1), Eq. (22) yields the asymptotic behavior 

G( I /;^) = |r(H)(|)' 2M + 0(i) 

and 

G(1 + 2 u- z) = Ir(l + 2 i /) + O (2 2 < 1+2 ")) , 

as z — ► 0. From Eq. (33) it follows that T int (L) vanishes 
for L -> 0 as 


'nnt(F) 


tq r(i + v) / l \ 
4 (l + 2*/)r(|i/|) \2Lo) 


2H 


(A4) 


In this limit the product T mt (L)a 2 (L) approaches a con- 
stant : 


Jim i r mt (L)a 2 (L) ~ 

Li — ►O 


7 or 0 r(l + v) 

2 1 + 2v * 


(A5) 


Appendix B: Effect of Finite Cut-off in 
the Evaluation of the Integral 
Correlation Time 

In this appendix we obtain an expression for the 
fractional error in the integral correlation time Tj n t pre- 
dicted by the spectral model for a circular disk of radius 
a as a function of the lag r max at which numerical in- 
tegration of the correlation function is stopped. The 


derivation is a simple extension of the calculations de- 
scribed in Bell and Kundu [2002]. 

The covariance of instantaneous area-averaged rain 
rate averaged over a circular disk of radius a separated 
by a time interval r is 


C aa (Q,t) 


(R' a (t)R' a ( 0)> 

-jj jf jf d 2 y <*(*, t)R'{ y, 0)) 

i// x // yc(x ~ y,r) 


J d 2 k J dwe 1 * ^ y ^c(k,o;) . 


Since c(k, a;) does not depend on the direction of k, 
the areal integrals can be done using the identity 


ra /*Z7T 

/ d 2 xe lkoc = I rdr l d<j>e 
Ja Jo Jo 

= 27ra 2 Ji(ka)/(ka) , 


ikr cos 4> 


where J\{x) is the Bessel function of the first kind of 
order 1. The integral over w with c(k,u;) as in Eq. (9) 
yields 


/ oo 

du>c{k,w) = V2nc(k,0) = nFoTk 

-OO 


After some algebra, one obtains 


Caa(0,t) — 

47or(l + i/) f°°dK Ji(kql) 
a 2 J 0 k v(k) 

x e -l r / r o|t/(K) _ 

(Bl) 

with the definitions 



a 

= a/L 0 , 

(B2) 

and 

»(*) 

= (1 + k 2 ) 1+v . 

(B3) 

For a circular disk of area L 2 the radius is a = Ljyfn* 


Using (Bl) one can easily show that 


^int 


f™drC AA (0,T) 

" f°°d£ 

JO k J 2 (k;) _ 
T ° f°°dK 


(B4) 
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Fractional error as function of cut-off 



ary, linear, first order autoregressive [AR(1)] process, 
evaluated from a single time series of finite length T. 

Consider the time series (x(t)} of length T generated 
by a single realization of a random variable X(t) under- 
going a linear stationary AR(1) process for which the 
lagged covariance function 

7(r) = <X'(i + T)X'(t)> (Cl) 

has the simple exponential form 

7(t) = a 2 e~W /TA . (C2) 

Its estimate, obtained from a finite time series, namely, 

c(r) = fJ o dtx'(t + T)x'(t ) 
can be thought of as a random variable 


Figure Bl. Plot of the fractional error as function 
of the ratio r max /ro computed for the typical spectral 
model parameters v = —0.25 and Lo = 70 km and three 
values of disk radius a. The curves marked (1), (2) and 
(3) represent the results for the radii a — 1 km, 10 km 
and 100 km. The corresponding true integral correla- 
tion times are 0.052r o , 0.19 t o and 0.65ro respectively. 


If we estimate r\ nt by introducing a finite upper cut-off 
in the r-integral, then the fractional error e(ro;T max ) 
defined as 

w) = - n - t ~ - - - n - . (B5) 

Tint 

is given by the formula 


e ( T 0 > 7 m ax) 


■C ex P [-( w/n>M«)] 

r°°dK 

Jo K V 2 (k) 


(B6) 

Eq. (B6) is numerically evaluated using Mathematica 
[Wolfram, 1999] to estimate the correction to the inte- 
gral correlation time discussed in section 3.3. A plot of 
e against the ratio r max /ro is shown in Figure Bl for 
typical parameter values. 


Y(T)~±£dtX'(t + T)X'(t) (C3) 

depending on the realization of X(t). We are interested 
in estimating the random error for the integral correla- 
tion time represented by the random variable 

r'T'max 

T in t = / Z(r)dr (C4) 

Jo 

where the variable Z(r) = y(r)/F(0) represents the 
lagged correlation function. The desired error estimate 
is simply the square root of the variance 

O ' 2 [Tint] = ((Tint - (Tint)) 2 ) , 

= (T 2 nt> - (Tint) 2 - 

The variance a 2 [Tint] can ^ en be expressed in the 
form 

/*T ma x / , %ax 

a 2 [T int ] = / / dT X dT^Z\T X )Z\T 2 )) . (C5) 

Jo Jo 

For simplicity of notation we write 

Y(t) = A,Y(0) = B. 


Appendix C: Estimate of Random Error 
for the Integral Correlation Time of an 
AR(1) Process 

In this appendix we obtain a simple estimate of the 
random error of the integral correlation time of a nor- 
mally distributed random variable undergoing a station- 


These random variables can each be expressed as the 
sum of respective mean values and fluctuations: 

A — a-j-Q:, J3 = 6 + /3, 

where 


a = (A) = 7 (t) , b = (B) = <x 2 . 
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Since up to second order in the fluctuations 


Eq. (C7) then reduces to 




it follows that 


where the cubic and higher order terms are neglected. 
In this approximation 


(Z'WM) = (+•#)-(+)(#) 

1 / I3 2 \ 

~ p- ( 0 : 10:2 + , (C6) 

where the subscripts 1 and 2 indicate evaluation at lags 
Ti and t 2 respectively. 

Eq. (C6) can be recast in a much simpler form if 
X(t) (and therefore X f (t)) is a gaussian random vari- 
able. Using the relation [see, e.g., Anderson , 1958] 


1 f°° 

(a xOt 2 ) dr ' h( r ' + T1- r 2 h(r / ) 

4* 7 (t ; + ri) 7 (r / - r 2 )]. (C9) 

Specializing Eq. (C9) to the case n = r 2 = 0, we also 
have 

9 r°° 

</? 2 > « ^ / <*r' 7 V). (CIO) 

Making use of the explicit form of the lagged covariance 
function 7(7*) given by Eq. (C2) and carrying out the 
remaining integration, we get 


(«1« 2 ) » 

Y [(In -r 2 | + r A )e |Tl ' 

~T2\/t A 


+ (n + r 2 + 7vO e -< Tl+T2 >M , (Cll) 

</J 2 ) a 

2<t 4 t a 

T ' 

(C12) 


Eq. (C5) then becomes 


{X\t x )X\t 2 )X\tz)X\U)) = 

(X\t x )X\ti))(X\t z )X\U)) 

+ {X'{t^X\t z ))(X'{t 2 )X\U)) 

+ {X\t{]X\U))[X\t 2 )X\ti)) , 

expressing the fourth moment of such a variable in 
terms of products of the second moment, the quantity 


a 2 [T int ] 


I 

T 



T*max 

dridr 2 


X [(In-ral + r^e-^-^x 
+ (n + t 2 + r A )e~ ( - Tl+T ^ /rA j 



/*T ma x 

/ d 
Uo 


n 2 


^ r ' e l T l/ r A 


(C13) 


<ai« 2 ) = (r'fajr'fa)) 

= (y(n)y(r 2 ))-(y(r 1 )>(y(r 2 )> 

can be reduced to the form 

1 f T f T 

( Q i Q 2> = jv J 0 J 0 dtldt 2 “ * 2 + n ~ 7 ’ 2 )^ 1 _ 

+ 7(<i -t 2 + Tib(f 1 - <2 - r 2 )l . (( 


The double integrals in the first term can be simplified 
using Eq. (C8) with r max instead of T . In the approxi- 
mation r max » To, we finally obtain the desired formula 

o' 2 [Tint] * [l + o (r A /r max )] , (C14) 

) T 

which is Eq. (42), since To = ta for the exponentially 
decaying correlation (C2). 


The double integrals can be easily performed with the 
help of the integral identity 

dtldt 2 m - 1 2) = J T dr (T - |t|)/(t) , (C8) 

which, for a function f(r) decaying exponentially as 
r — > 00, in the approximation T » ta further simplifies 
to 

dt\dt 2 f(ti — f 2 ) « T f dr' f{r') . 

J — OO 
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within about 1% of what was reported in BK96 from a 
direct fit to the lagged autocorrelation function. 


4. Results and Discussion 

The full set of parameter values for which the model 
best fits the TOGA CO ARE data are given in Table 1, 
along with the corresponding values for GATE Phase I 
determined in BK96 using a somewhat less exact (and 
more arduous) trial-and-error curve fitting method in- 
stead of the analytical method adopted here. 

4.1. Model Parameters v , Lo and 70 

Table 1 summarizes the TOGA COARE model pa- 
rameters, We see that the values of the power law expo- 
nent v range between about —0.21 to -0.34 indicating 
a much stronger singular behavior than was found for 
the GATE Phase I data (v = —0.11). The character- 
istic time scales tq lie in the range of roughly 4.5 to 
8.2 hours, considerably shorter than the GATE Phase 
I value ro « 13 hours. The characteristic length Lo 
ranges between about 54 and 94 km indicating a some- 
what sharper fall-off of the spatial correlations com- 
pared to GATE Phase I for which Lo = 104 km. The 
results exhibiting the quality of the fit are shown in Fig- 
ure 2 for the TOGA and MIT radars. It is seen that in 
most cases the fit is quite good up to grid size L — 64 
km, beyond which the approximation on which Eq. (24) 
is based is expected to break down. Our results seem to 
be consistent with the model prediction that the vari- 
ance <J 2 {L) has a power law behavior <r 2 (L) ~ L~ as 
L — > 0. Robustness of the parameterization was checked 
by comparing the variance of area-averaged rain rate 
(<x 2 )th for an area A — L 2 with L — 128 km predicted 
by the exact formula (20), rather than the asymptotic 
formula (24), against the observed variance (<7^) 0 bs for 
the 128 km box centered at the two ship locations. Ta- 
ble 2 shows the comparisons, which are quite good. 

For GATE data we find that for L between 4 km and 
56 km, the variance cr 2 {L) is quite accurately fitted (not 
shown here) by the formulas 

2 / r \ f 17.3L~ 0 24 — 4.66 mm 2 h“ 2 Phase I 

a 1 ' \ 15.0L- 0 38 - 1.92 mm 2 h -2 Phase II 1 

which through Eqs. (24) to (26) yield the parameter 
values 

7o = 1-02 mm 2 h" 2 , v — -0.12, Lo = 93.8 km 


Variance as function of area for TOGA COARE 




Figure 2. Comparison between the observed variance 
of area-averaged rain rate a 2 (L) for the TOGA COARE 
data and the corresponding spectral model predictions 
for various box sizes L ranging from 2 km to 128 km 
using Eqs. (24)-(26) and the parameter values listed in 
Table 1. 


for Phase I and 

7o — 0.63 mm 2 h“ 2 , v = -0.19, L 0 = 82.7 km 




